R Notebook

What are the most frequent minor variants in the high-confidence samples?

source("./scripts/startup.R")

The following packages are a base install and will not be unloaded:



The following packages were not previously loaded:



Loading required package: pacman

patient_var_30 = read_feather("processing/patient_var_30.arrow")
maf_histogram = patient_var_30 %>% ggplot(aes(minorfreq)) + geom_histogram(binwidth = 0.01) + theme_pubr()
quantile(patient_var_30$minorfreq, probs = c(.25,.5,.75), na.rm = T)

       25%        50%        75% 
0.01252470 0.01717647 0.02788318 

ggsave("ggsave/maf_histogram.pdf", maf_histogram, width = 3, height = 3)

patient_var_30_for_rank = patient_var_30 %>% drop_na(gene) %>% filter(gene !="") %>%
  mutate(label = paste0(ref_sym, aapos), 
         refnt = str_sub(ref_codon, start = codon_pos, end = codon_pos))
bp = c("A", "T", "C", "G")
lineage_defining_mutations = fread("20220103-TRACE-LineageDefinitions-v9.1.txt", data.table = F) %>%
  filter(nt_ref %in% bp & 
           nt_alt %in% bp & 
           variation_type == "SNP") %>% select(nt_pos, nt_ref, nt_alt)
ldm_map = patient_var_30_for_rank %>% 
  mutate(ldm  = ifelse(ntpos %in% lineage_defining_mutations$nt_pos, TRUE, FALSE)) %>% 
  select(gene, label, ldm) %>% distinct %>% group_by(gene, label) %>% arrange(-ldm) %>%
  filter(row_number()==1)
  

highly_shared_sites_ranked = patient_var_30_for_rank %>% drop_na(gene) %>% filter(gene !="") %>%
  mutate(label = paste0(ref_sym, aapos), 
         refnt = str_sub(ref_codon, start = codon_pos, end = codon_pos)) %>%
  group_by(mcov_id, gene, label) %>%
  tally() %>% 
  arrange(-n) %>% ungroup %>% group_by(gene, label) %>% 
  summarize(label_count = n()) %>% ungroup() %>% arrange(-label_count) %>% mutate(rank = 1:nrow(.)) %>%
  left_join(ldm_map) %>% unique()

# highly_shared_sites_ranked_50 = highly_shared_sites_ranked %>% filter(rank <= 25)
# 
# ((
#   hss_maf_50<-patient_var_30 %>% filter(paste0(gene, aapos) %in% paste0(highly_shared_sites_ranked_50$gene, parse_number(highly_shared_sites_ranked_50$label))) %>% 
#     ggplot(aes(x = as.factor(ntpos), y = minorfreq)) + 
#     #geom_point(alpha = 0.2, position = position_jitter(width = 0.1)) + theme_bw() + 
#     geom_boxplot(outlier.shape = NA, alpha =0.5) + theme_bw() +
#     theme(axis.text.x = element_text(angle=90)) + #geom_boxplot(outlier.alpha = 0) + 
#     xlab("Nucleotide position") + ylab("Minor allele freq") ))
# 
# ggsave("ggsave/hss_maf_50.pdf", hss_maf_50, width = 5, height = 3)


#%>% 
 # left_join(patient_var_30 %>% select(ntpos, gene, ref_sym, aapos, ref_codon, codon_pos) %>% distinct)

# ((
#   highly_shared_sites<- patient_var_30 %>% group_by(ntpos) %>% tally() %>% arrange(desc(n)) %>% 
#        drop_na(ntpos) %>% slice_max(n, n = 35) %>% pull(ntpos) %>% as.numeric()
# ))


sample_n = length(unique(patient_var_30$mcov_id))
#for (this_gene in gene_limits$gene_id) {
  plot_ranks = highly_shared_sites_ranked %>% filter(rank <= 35) %>% 
    ggplot(aes(rank, label_count, label = label)) + 
  geom_line(data = highly_shared_sites_ranked %>% 
              select(!gene) %>% filter(rank <= 35), aes(rank, label_count), color = "grey") +
  scale_y_continuous(trans = "log2", breaks = 2^seq(0,12,1),
                     sec.axis = sec_axis(~./sample_n, labels = scales::label_percent(),
                                           breaks = 2^seq(0,8,1)/100)) + theme_pubr() +
  geom_point(aes(rank, label_count), color = "grey") +
  geom_text_repel(aes(label = label, color = ldm), ylim = c(6,11), xlim = c(NA, NA), angle = 90,
                  segment.size = 0.6, segment.curvature = -1e-20, segment.linetype = 3,
                  max.overlaps = Inf) + coord_cartesian(clip = "off") +
    scale_color_manual(values = c("salmon", "grey")) +
    facet_wrap(~gene)
  
  plot_ranks
#}
ggsave("ggsave/plot_ranks.pdf", plot_ranks, height = 5.5, width = 7)

# Across genome


lineage_defining_mutations

plot_data_LDM_prop = patient_var_30 %>% group_by(MCoVNumber) %>% drop_na(ntpos) %>% 
  summarize(iSNV_count = n(), iSNV = paste0(ntpos, collapse = ","), 
            iSNV_LDM_count = length(ntpos[ntpos %in% lineage_defining_mutations$nt_pos]),
            iSNV_LDM = paste0(ntpos[ntpos %in% lineage_defining_mutations$nt_pos],
                              collapse = ","), prop = iSNV_LDM_count/iSNV_count)

plot_data_LDM_prop %>% ggplot(aes(as.factor(iSNV_count), prop)) + geom_violin(draw_quantiles=c(0.5))

patient_var_30_reversions = patient_var_30 %>% mutate(reversion = minor == str_sub(ref_codon, 
                            start = codon_pos, end = codon_pos)) %>% drop_na(reversion)

# How many of our alleles were reversions?
patient_var_30_reversions %>% summarize(proportion_reversion = sum(reversion==T)/
                                          nrow(patient_var_30))

# How many of our alleles that are reversions are at LDM sites?
patient_var_30_reversions %>% mutate(ldm = ifelse(ntpos %in% lineage_defining_mutations$nt_pos,1,0)) %>% 
  summarize(proportion_reversion = sum(ldm == 1 & reversion == T)/sum(reversion==T))

# obscure the consensus in the VCF
problem_sites_global = fread("problematic_sites_sarsCov2_v8-20211027.vcf", skip = 88) %>%
  filter(FILTER != "mask") %>%
  filter(!grepl('single_src|nanopore', INFO)) %>% pull(POS)
problem_sites_houston = fread("problematic_sites_sarsCov2_v8-20211027.vcf", skip = 88) %>%
  filter(grepl('Houston', INFO)) %>% pull(POS)
problem_sites = unique(c(problem_sites_global, problem_sites_houston))
problem_sites

 [1]    76    78   320   538   660  1001  1814  1947  2087  5393  5498  6309  6310  6312  6866  7396
[17] 12685 13512 13513 13514 13686 13693 18505 18506 19338 19344 19369 19732 21302 21304 21305 21658
[33] 22393 22410 22488 22515 23144 24673 25798 26709 27792 28881 28882 28883 29378  8658 12698 15103
[49] 16130 16132 22892

# UNCOMMENT BELOW IF NECESSARY TO REGENERATE DATA
# consensus <-fread("consensus_minor_changes_20220713.csv", data.table = F) %>% filter(!ntpos %in% problem_sites) %>% filter(!ntpos %in% c(1:265)) %>% filter(!ntpos>29674)
# 
# 
# 
# consensus 
# 
# #
# consensus_count = consensus %>% 
#   mutate(MCoVNumber = mcov_reformat(name)) %>% 
#   select(MCoVNumber, major, ntpos, refnt) %>% 
#   filter(major != refnt) %>% 
#   unique() %>%
#   group_by(ntpos) %>%
#   summarize(MCoVNumber = unique(MCoVNumber)) 
# 
# runs = fread("sample_date_and_run.csv", data.table = F) %>%
#     mutate(MCoVNumber = mcov_reformat(mcov_id))
# 
# consensus_runs = consensus_count %>% left_join(runs) %>% drop_na() %>%
#   select(ntpos, MCoVNumber_possible_contaminant = MCoVNumber, run_consensus = run_group) #has each of the consensus mutations with the MCOV
# 
# # uncomment below if necessary
# #patient_var_30_contaminated = patient_var_30 %>% left_join(consensus_runs) # intensive
# #write_feather(patient_var_30_contaminated, "processing/patient_var_30_contaminated.arrow")
# 
patient_var_30_contaminated = read_feather("processing/patient_var_30_contaminated.arrow")
contaminated_df = patient_var_30_contaminated %>% mutate(ref_nt = str_sub(ref_codon,
                            start = codon_pos, end = codon_pos)) %>%
  select(MCoVNumber, run_group, ntpos, ref_nt, major, minor, MCoVNumber_possible_contaminant, run_consensus)

contaminated_df_tallied = contaminated_df %>%
  filter(run_group == run_consensus) %>%
  group_by(MCoVNumber) %>%
  summarize(n_var_contam = length(unique(ntpos)), MCoVNumber_possible_contaminant) %>%
  group_by(MCoVNumber, MCoVNumber_possible_contaminant) %>%
  summarize(n_var_contam, single_sample = n(), prop_single_sample = single_sample / n_var_contam)

`summarise()` has grouped output by 'MCoVNumber'. You can override using the `.groups` argument.`summarise()` has grouped output by 'MCoVNumber', 'MCoVNumber_possible_contaminant'. You can override using the `.groups` argument.

contaminated_df_tallied_top = contaminated_df_tallied %>% group_by(MCoVNumber) %>% top_n(n=1) %>%
  filter(row_number()==1)

Selecting by prop_single_sample

#Note prop_single_sample is single_sample / n_var_contam (proportion of number of minor variants that resememble consensus in another sample / number of minor variants that resemble any consensus in the run)
# prop single contam is the same as above but the denom is out of n_var
patient_counts_30 = read_feather("processing/patient_counts_30.arrow")
patient_counts_30_anno = patient_counts_30 %>% left_join(contaminated_df_tallied_top) %>%
  mutate(run_num = as.numeric(str_replace(run_group, "Run_", ""))) %>%
  replace(is.na(.),0)

Joining, by = "MCoVNumber"

# Plot: n_var x prob that it's from a single source
plot_contam = patient_counts_30_anno %>% mutate(prop_contam = n_var_contam/n_var, prop_single_contam = single_sample/n_var) 

plot_contam %>% select(n_var, admitted_hospital, prop_contam, prop_single_contam) %>% 
  pivot_longer(cols= c(prop_contam, prop_single_contam), names_to = "type", values_to = "prop")  %>%
  ggplot(aes(x = admitted_hospital, y = prop, 
             fill = type)) + 
  geom_violin(draw_quantiles = c(0.25,0.5,0.75))

median_contam_contam = quantile(plot_contam %>% filter(n_var > 0) %>% 
                                  pull(prop_contam), probs = c(0.5))
print(median_contam_contam)

50% 
0.5 

plot_contam_contam = plot_contam %>% 
  select(n_var, admitted_hospital, prop_contam, prop_single_contam) %>% 
  ggplot(aes(x = n_var, y = prop_contam, 
             fill = admitted_hospital, color = admitted_hospital)) + 
  geom_point(alpha = 0.02) + geom_smooth() + theme_pubr() +
  geom_hline(yintercept = median_contam_contam, linetype = "dashed")

plot_contam_contam = ggMarginal(plot_contam_contam, groupColour = T, groupFill = T, type = "violin",
           draw_quantiles = c(0.25, 0.5, 0.75))

`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'Warning: Removed 450 rows containing non-finite values (`stat_smooth()`).`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'Warning: Removed 450 rows containing non-finite values (`stat_smooth()`).Warning: Removed 450 rows containing missing values (`geom_point()`).

median_single_contam = quantile(plot_contam %>% filter(n_var > 0) %>% 
                                  pull(prop_single_contam), probs = c(0.5))
print(median_single_contam)

      50% 
0.3333333 

plot_contam_single_contam = plot_contam %>% 
  select(n_var, admitted_hospital, prop_contam, prop_single_contam) %>% 
  ggplot(aes(x = n_var, y = prop_single_contam, 
             fill = admitted_hospital, color = admitted_hospital)) + 
  geom_point(alpha = 0.02) + geom_smooth() + theme_pubr() + 
  geom_hline(yintercept = median_single_contam, linetype = "dashed")

plot_contam_single_contam = ggMarginal(plot_contam_single_contam, 
                                       groupColour = T, groupFill = T, 
                                       type = "violin",
           draw_quantiles = c(0.25, 0.5, 0.75))

`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'Warning: Removed 450 rows containing non-finite values (`stat_smooth()`).`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'Warning: Removed 450 rows containing non-finite values (`stat_smooth()`).Warning: Removed 450 rows containing missing values (`geom_point()`).

plot_contam_arranged = ggarrange(plot_contam_contam, plot_contam_single_contam,
                                 align = "v",
                                 labels = "AUTO")
plot_contam_arranged
ggsave("ggsave/plot_contam_arranged.pdf", plot_contam_arranged, height = 3, width = 6)

#consensus_sites_lowct %>% filter(refnt == major)
patient_counts_30 = read_feather("processing/patient_counts_30.arrow")


# criteria to count the mutation
# consensus = if ntpos matches in the list above, then 
# the minor will also match the nt_ref (i.e. reversion) or the nt_alt


ldm = c(paste0(lineage_defining_mutations$nt_alt,
                                      lineage_defining_mutations$nt_pos,
                                      lineage_defining_mutations$nt_ref),
                               paste0(lineage_defining_mutations$nt_ref,
                                      lineage_defining_mutations$nt_pos,
                                      lineage_defining_mutations$nt_alt))
        
minor_consensus_plotdata = patient_var_30 %>% #left_join(consensus) %>% 
    mutate(ref_nt = str_sub(ref_codon, 
                            start = codon_pos, end = codon_pos),
      ref_mutation = paste0(ref_nt, ntpos, minor),
      mutation = paste0(major, ntpos, minor), consensus = 
             ifelse(ntpos %in% lineage_defining_mutations$nt_pos, 
                    ifelse((mutation %in% ldm) | (ref_mutation %in% ldm),
                           TRUE, FALSE),      
                           FALSE)) %>% #filter(consensus == FALSE) %>%
  group_by(ntpos, consensus) %>% 
  summarize(consensus = consensus, n = n()) %>%
  unique %>% drop_na %>% ungroup %>% arrange(-n) %>%
  mutate(rank = 1:nrow(.))

`summarise()` has grouped output by 'ntpos', 'consensus'. You can override using the `.groups` argument.

minor_prevalence_across_genome_unlog <- minor_consensus_plotdata  %>%
    ggplot(aes(x = ntpos, y = n, fill = consensus, color = consensus)) + 
    # label your data too for the SALMON
    geom_bar(stat = "identity") + 
  scale_color_manual(values = c("salmon","grey")) + scale_fill_manual(values = c("salmon","grey")) +
    
      geom_point( aes(x = ntpos, y= n), shape = "") +
    #scale_color_manual(values = c("salmon","black")) + 
  theme_bw() + 
  geom_hline(yintercept = nrow(patient_counts_30)*0.01, linetype="dotted") + 
  xlab("Nucleotide position") + ylab("No. samples w/minor variant at site") + 
  annotate("rect", xmin=27894, xmax=28295, ymin=0, ymax=Inf, alpha=0.2, fill="#85D4E3") + 
  annotate("rect", xmin=28274, xmax=29533, ymin=0, ymax=Inf, alpha=0.2, fill="#F4B5BD") + 
  annotate("rect", xmin=21563, xmax=25384, ymin=0, ymax=Inf, alpha=0.2, fill="#FAD77B") +  xlim(0,29903)

minor_prevalence_across_genome_unlog

minor_prevalence_across_genome_unlog_labeled = 
  minor_prevalence_across_genome_unlog + 
  geom_text_repel(max.overlaps = Inf, ylim  = c(500,NA), size = 3, angle = 90, 
                  segment.linetype = 3, segment.size = 0.6, 
                  segment.curvature = -1e-20, 
                  arrow = arrow(length = unit(0.015, "npc"), 
                                type = "closed"),
                  data = minor_consensus_plotdata %>%
                    filter(rank < 50),
                    aes(x = ntpos, y = n, label = ntpos, color = consensus, fill = consensus))

Warning: Ignoring unknown aesthetics: fill

minor_prevalence_across_genome_unlog_labeled

genes <- fread("ntpos_gene_update.csv", data.table = F)

gene_limits<- genes %>% group_by(gene_id) %>% summarise(start=min(ntpos), end=max(ntpos)) %>% filter(gene_id!="") %>% mutate(molecule="")

gene_arrows<-ggplot(gene_limits, aes(xmin = start, xmax = end, y=molecule, label = gene_id)) +
  geom_gene_arrow(arrowhead_height = unit(3, "mm"), arrowhead_width = unit(1, "mm")) + geom_gene_label() + geom_segment(aes(x=266,y=1.5,xend=13468,yend=1.5), size=0.2) + xlim(0,29903) +
  annotate("text", label="ORF1a", x=5000, y=1.41, size=3) + geom_segment(aes(x=13468,y=1.4,xend=21555,yend=1.4), size=0.2) + annotate("text", label="ORF1b", x=18000, y=1.31, size=3) +
  theme_genes() + ylab(NULL) + xlab(NULL) 

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.

((
  genome_fig_unlog<-plot_grid(gene_arrows, minor_prevalence_across_genome_unlog_labeled + theme(legend.position = "bottom"), nrow=2, rel_heights=c(0.25,1), axis="lr", align="hv")
))



ggsave("ggsave/genome_fig_unlog.pdf", plot = genome_fig_unlog, height = 5, width = 12)

Characteristics of highly recurrent minor variants

# What's the range of minor allele frequencies the highly recurrent minor variants are found at?
highly_shared_sites_top_ranked = highly_shared_sites_ranked %>% filter(rank <= 35)
highly_shared_sites = patient_var_30_for_rank %>% filter(paste0(gene,".",label) %in%
                                     paste0(highly_shared_sites_top_ranked$gene, ".", 
                                            highly_shared_sites_top_ranked$label)) %>%
          group_by(label, ntpos) %>% 
  summarize(ntpos, ntpos_count = n()) %>% distinct() %>% 
  arrange(-ntpos_count) %>% 
  filter(ntpos_count > (patient_var_30$MCoVNumber %>% unique %>% length)*0.01) %>% pull(ntpos)

`summarise()` has grouped output by 'label', 'ntpos'. You can override using the `.groups` argument.

((
  hss_maf<-patient_var_30 %>% filter(ntpos %in% highly_shared_sites) %>% 
    ggplot(aes(x = as.factor(ntpos), y = minorfreq)) + 
    geom_point(alpha = 0.2, position = position_jitter(width = 0.1)) + theme_bw() + 
    geom_boxplot(outlier.shape = NA, alpha =0.5) + 
    theme(axis.text.x = element_text(angle=90)) + #geom_boxplot(outlier.alpha = 0) + 
    xlab("Nucleotide position") + ylab("Minor allele freq") + 
    scale_y_continuous(trans="log2", breaks = c(0.01*2^(0:5), 0.5))
  ))

ggsave("ggsave/hss_maf.pdf", width = 4, height = 1, plot = hss_maf)

((
  hss_maf_unlog<-patient_var_30 %>% filter(ntpos %in% highly_shared_sites) %>% 
    ggplot(aes(x = as.factor(ntpos), y = minorfreq)) + 
    geom_point(alpha = 0.2, position = position_jitter(width = 0.1)) + theme_bw() + 
    geom_boxplot(outlier.shape = NA, alpha = 0.5) + 
    theme(axis.text.x = element_text(angle=90)) + #geom_boxplot(outlier.alpha = 0) + 
    xlab("Nucleotide position") + ylab("Minor allele freq") #+ 
    #scale_y_continuous(trans="log2", breaks = c(0.01*2^(0:5), 0.5))
  ))
ggsave("ggsave/hss_maf_unlog.pdf", width = 4, height = 2, plot = hss_maf_unlog)

#How many different changes do you find at each recurrently-mutated position? How many such mutations are a reversion to the reference?
hss_info <- fread('hss_data2.csv', data.table = F)

((
  hss_plot = patient_var_30_for_rank %>% mutate(label_gene = paste0(gene, ": ", label)) %>%
    filter(ntpos %in% highly_shared_sites) %>% 
    group_by(ntpos, major, minor) %>% mutate(n=n(), median_frequency = median(minorfreq)) %>%     mutate(change=paste0(major,">",minor)) %>% 
    mutate(ref_reversion=if_else(minor==refnt,"yes","no")) %>% 
    mutate(median_MAF=if_else(median_frequency<0.02, 
                              "low (median <2% MAF)","high (median >=2% MAF)")) %>%
    select(change, n, ref_reversion, ntpos, label_gene, median_MAF) %>% distinct %>%
    ggplot(aes(x=change, y=n, color=ref_reversion)) + 
    geom_bar(stat="identity", aes(fill=median_MAF)) + 
    scale_color_manual(values=c("white","red")) + 
    scale_fill_manual(values=c("darkmagenta","steelblue")) + 
    facet_wrap(ntpos~label_gene, scales="free") + 
    theme_pubr() + theme(axis.text.x=element_text(angle=90))
))

Adding missing grouping variables: `major`, `minor`

ggsave('ggsave/hss_plot.pdf', plot = hss_plot, width = 11, height = 13)

ldm_vector = lineage_defining_mutations$nt_pos
heatmap_spectra_tmp = patient_var_30_for_rank %>% 
  mutate(label_gene = paste0(gene, ": ", label)) %>%
    filter(ntpos %in% highly_shared_sites) %>% 
    group_by(ntpos, major, minor) %>% mutate(n=n(), median_frequency = median(minorfreq)) %>%
  mutate(change=paste0(major,">",minor), 
         transition = ifelse(change %in% c("C>T","T>C","A>G","G>A"), 1, 0)) %>%
     mutate(ref_reversion=if_else(minor==refnt,1,0)) %>% ungroup() %>%
  select(change, n, ntpos, label_gene) %>% distinct %>% 
  group_by(ntpos) %>% mutate(ntpos_n = sum(n)) %>% ungroup() %>%
  mutate(n_prop = n/ntpos_n) %>% select(-n) %>% 
  mutate(ldm = as.factor(ifelse(ntpos %in% ldm_vector, 1, 0))) %>%
  mutate(label = paste0(label_gene, " - ", ntpos, " (", ntpos_n, ")"))


annotation_row_df = heatmap_spectra_tmp %>% select(label, ldm) %>% 
  distinct %>% as.data.frame() %>%
  column_to_rownames("label")
annotation_row_df

#%>%
heatmap_spectra = heatmap_spectra_tmp %>% select(-ldm) %>% spread(change, n_prop) %>% 
  arrange(-ntpos_n) %>% distinct() %>%  
  select(-c(ntpos, label_gene, ntpos_n)) %>% replace(.,is.na(.),0) %>%
  column_to_rownames("label")

annotation_col_df = heatmap_spectra_tmp %>% select(change) %>% distinct %>% 
  mutate(transversion = 
           as.factor(ifelse(!(change %in% c("C>T","T>C","A>G","G>A")), 
                            1, 0))) %>% 
  column_to_rownames("change")

maf = heatmap_spectra_tmp %>% select(ntpos, label) %>% distinct %>%
  left_join(patient_var_30 %>% select(ntpos, minorfreq)) %>% 
  select(label, minorfreq)

Joining, by = "ntpos"

maf_list = split(maf$minorfreq,maf$label)

row_ha = rowAnnotation(df = data.frame(
  ldm = annotation_row_df[rownames(heatmap_spectra),]), 
  col = list(ldm = c(`1`="grey", `0`="white")),
  MAF = anno_boxplot(maf_list, pch = 20, size = unit(0.5, "mm")))


heatmap_spectra_reorder = heatmap_spectra %>% select(c("A>G","G>A","C>T","T>C",
                                               "A>C", "C>A", "A>T", "T>A",
                                               "C>G", "G>C", "G>T", "T>G"))
col_ha = columnAnnotation(df = data.frame(
  transversion = annotation_col_df[colnames(heatmap_spectra_reorder),]), 
                          col = list(transversion = c(`1`="grey", `0`="white")))

heatmap_spectra_hss = Heatmap(heatmap_spectra_reorder,
                              name = "prop/site",
                              col = colorRampPalette(
                                c("white", "brown4"))(100),
                              show_row_dend = F, 
                              cluster_columns = F, left_annotation = row_ha,
                              rect_gp = gpar(col = "grey90", lwd = 1),
                              top_annotation = col_ha,
                              heatmap_legend_param = 
                                list(direction = "vertical")) 

Warning: The input is a data frame-like object, convert it to a matrix.Warning: Values in row annotation 'MAF' have a different order of names from the matrix
row names. It may lead to wrong conclusion of your data. Please double check.

pdf(qq("ggsave/heatmap_spectra_hss.pdf"), width = 6.5, height = 7)
draw(heatmap_spectra_hss, merge_legend = TRUE, 
     heatmap_legend_side = "right", 
    annotation_legend_side = "right")
dev.off()

null device 
          1 

heatmap_spectra_hss

# 
#           border_color = "grey90", colorRampPalette(c("white", "brown4"))(100),
#          treeheight_row = 0, treeheight_col = 10, annotation_col = annotation_col_df,
#          annotation_row = annotation_row_df, annotation_colors = annotation_colors)
# pacman::p_load(ComplexHeatmap)
         

#all minor variants that are found in at least 2% of high-coverage samples in Lythgoe et al https://www.science.org/doi/suppl/10.1126/science.abg0821/suppl_file/abg0821-lythgoe-sm.pdf

Warning message:
The closing backticks on line 240 ("```") in 03_recurrent_mutations_analysis.Rmd do not match the opening backticks " ```" on line 209. You are recommended to fix either the opening or closing delimiter of the code chunk to use exactly the same numbers of backticks and same level of indentation (or blockquote). 

lythgoe_hss<-c(29320, 25628, 25807, 20364, 357, 25627, 25532, 11083, 239, 238, 356, 16740, 20989, 19393, 15009, 21826, 19937, 26289,28253,25507,29862,11052,29864,25529,369,22453,13571,29538,21575,75,20993,19473,15474,635,1226,29747,24933,19210,13303,28936,26334,16887)

tonkin_hss<-c(11075L, 11083L, 11074L, 522L, 26780L, 1547L, 14408L, 28253L, 13914L, 26137L, 635L, 241L, 26785L, 683L, 23403L, 28881L, 13778L, 25521L, 14805L, 9491L, 11071L, 13780L, 21575L, 29242L, 9430L, 203L, 3037L, 686L, 3096L, 16375L, 16466L, 17167L, 23559L, 9438L, 21137L, 29051L, 23555L, 26333L, 9474L, 13571L, 16887L, 26787L, 28603L, 24213L, 11651L, 27925L, 514L, 26781L, 558L, 1912L, 9479L, 11073L, 12578L, 26144L, 26681L, 29241L)
intersect(lythgoe_hss, highly_shared_sites)

[1] 11083 28253

intersect(tonkin_hss, highly_shared_sites)

[1] 11083 28253 14805

homoplasic_site_list<-c(187L, 1059L, 2094L, 3037L, 3130L, 6990L, 8022L, 10323L, 10741L, 11074L, 11083L, 13408L, 14786L, 15324L, 19684L, 20148L, 21137L, 21575L, 24034L, 24378L, 25563L, 26144L, 26461L, 26681L, 28077L, 28826L, 28854L, 29700L)
intersect(homoplasic_site_list, highly_shared_sites)

[1] 11083 25563

ggsave("ggsave/gisaid_vs_hmh_plot_marginal.pdf", plot = gisaid_vs_hmh_marginal, width = 6, height = 5)

Warning: ggrepel: 15 unlabeled data points (too many overlaps). Consider increasing max.overlaps

# MUTATION SPECTRA ANALYSIS

# reproducibility analysis from RMD 01
plot_nonreproducible_spectra_heatmap = readRDS(file = "plot_nonreproducible_spectra_heatmap.rds")
plot_reproducible_spectra_heatmap = readRDS(file = "plot_reproducible_spectra_heatmap.rds")

hmap<-table(patient_var_30$major[patient_var_30$major!=""], patient_var_30$minor[patient_var_30$minor!=""]) %>% data.frame() %>%
  ggplot(aes(x=Var1, y=Var2, fill=Freq/sum(Freq))) + geom_tile(colour = "black") +
  scale_fill_gradient(low = "white",
                      high = "darkblue") +
  theme_minimal() +   labs(fill = "Number",
       x = "Consensus allele",
       y = "Minor allele", title = "Post-sample filter (Ct, n_var)") #Most frequently C>T mutations

# Now beecause some peaks are really prevalent like the twin peaks at 29187/8, 
# it's important to not factor that in too much. So let's just assign those mutations
# to equate to the median number a mutation appears which is 1.

mutation_spectra = patient_var_30 %>% unite(col = "spectra_mutation", 
                                            major, ntpos, minor, remove = FALSE)
mutation_spectra_counts = mutation_spectra %>% count(spectra_mutation)

mutation_spectra_counts %>% ggplot(aes(y=n)) + geom_boxplot() + 
  scale_y_continuous(trans = "log1p")

mutation_spectra_counts$n %>% quantile() # median is 1

mutation_spectra_unique = mutation_spectra %>% 
  select(major, minor, spectra_mutation) %>% unique()

hmap_unique <-table(mutation_spectra_unique$major[mutation_spectra_unique$major!=""], mutation_spectra_unique$minor[mutation_spectra_unique$minor!=""]) %>% data.frame() %>%
  ggplot(aes(x=Var1, y=Var2, fill=Freq/sum(Freq))) + geom_tile(colour = "black") + # grid color
  scale_fill_gradient(limits = c(0,0.4), low = "white",
                      high = "darkblue") +
  theme_minimal() +   labs(fill = "Fraction",
       x = "Consensus allele",
       y = "Minor allele", title="Unique minor variants")

rescale_fill = function(hmap_unique) {
  hmap_unique_rescaled = hmap_unique + 
    scale_fill_gradient(limits = c(0,0.4), low = "white",
                      high = "darkblue") +
    geom_label(fill = "white", alpha = 0.3, aes(x = Var1, y = Var2, label = round(Freq/sum(Freq), digits = 2)))
  return(hmap_unique_rescaled)
}

((
  rep_nucleotides = ggarrange(rescale_fill(plot_nonreproducible_spectra_heatmap), 
                              rescale_fill(plot_reproducible_spectra_heatmap), 
                              rescale_fill(hmap), 
                              rescale_fill(hmap_unique), common.legend = T)
))
rep_nucleotides

ggsave("ggsave/heatmap_spectra_replicate_variants.pdf", rep_nucleotides, height = 6, width = 6)